# TODO: Add comment
# 
# Author: rogb
###############################################################################

armod <- function(x,maxAR,reduce=FALSE,include.mean=FALSE,usePCA=FALSE,pVar=0.99){
	#maxAR <- 12
	n <- length(x)
	usedSeq <- n - maxAR
	idxOrg <- 1:usedSeq
	idxExpl <- t(t(matrix(rep((idxOrg),maxAR),ncol=maxAR)) + 1:maxAR)
	xr <- rev(x)
	if(usePCA){
		data <- data.frame(matrix(xr[idxExpl],ncol=maxAR))
		names(data) <- c(paste("MA",1:maxAR,sep=""))
		
		pc <- princomp(data)
		cumVar <- cumsum(pc$sdev^2/sum(pc$sdev^2))
		nPC <- min(which(cumVar>=pVar))
		data <- data.frame(xr[idxOrg],pc$scores[,1:nPC])
		names(data) <- c("Orig",paste("Comp",1:nPC,sep=""))
		mod <- lm(Orig~.-1,data=data)	
		mod <- step(mod,trace = FALSE)
		maCoefs <- t(pc$loadings[,1:length(mod$coefficients)] %*% as.matrix(mod$coefficients))
		data <- data.frame(xr[idxOrg],matrix(xr[idxExpl],ncol=maxAR))
		names(data) <- c("Orig",paste("MA",1:maxAR,sep=""))
		mod <- lm(Orig~.-1,data=data)
		coefNames <- names(mod$coefficients)
		mod$coefficients[] <- maCoefs
	}
	else{
		data <- data.frame(xr[idxOrg],matrix(xr[idxExpl],ncol=maxAR))
		names(data) <- c("Orig",paste("MA",1:maxAR,sep=""))
	}
	
	if(include.mean==TRUE){
		mod <- lm(Orig~.,data=data)	
	}
	else{
		mod <- lm(Orig~.-1,data=data,weights=w)		
	}
	if(reduce){
		mod <- step(mod,trace = FALSE)
	}
	mod
}

predict.armod <- function(mod){
	maUsed <- names(mod$coef)
	x <- mod$fitted.values + mod$residuals
	if(is.element("(Intercept)",maUsed)){
		if(length(maUsed)>1){
			maIdx <- as.numeric(substring(maUsed[2:length(maUsed)],3,nchar(maUsed[2:length(maUsed)])))
			newdata <- eval(parse(text=paste("data.frame(",paste(paste(maUsed[2:length(maUsed)],"=",x[maIdx]),collapse=","),")",sep="")))
			out <- predict(mod,newdata=newdata,se.fit =TRUE,interval = "prediction")
		}else{
			out <- predict(mod,newdata=data.frame(),se.fit =TRUE,interval = "prediction")
		}
	}else{
		maIdx <- as.numeric(substring(maUsed,3,nchar(maUsed)))
		newdata <- eval(parse(text=paste("data.frame(",paste(paste(maUsed,"=",x[maIdx]),collapse=","),")",sep="")))
		out <- predict(mod,newdata=newdata,se.fit =TRUE,interval = "prediction")
	}
	out
}
